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Abstract. This paper is concerned with fast and accurate computation of exterior wave 
equations truncated via exact circular or spherical nonreflecting boundary conditions (NRBCs, 
which arc known to be nonlocal in both time and space). We first derive analytic expressions 
for the underlying convolution kernels, which allow for a rapid and accurate evaluation of the 
convolution with O(Nt) operations over Nt successive time steps. To handle the nonlocality in 
space, we introduce the notion of boundary perturbation, which enables us to handle general 
bounded scatters by solving a sequence of wave equations in a regular domain. We propose 
an efficient spectral-Galerkin solver with Newmark's time integration for the truncated wave 
equation in the regular domain. We also provide ample numerical results to show high-order 
accuracy of NRBCs and efficiency of the proposed scheme. 

1. Introduction 

Wave propagation and scattering problems in unbounded media arise from diverse applica- 
tion areas such as acoustics, aerodynamics, electromagnetics, antenna design, oceanography and 
among others (see, e.g., [THl [Ml EH] ) ■ Various approaches have been proposed for their numer- 
ical studies that include the boundary element methods (cf. [7]), infinite element methods (cf. 
[13]), perfectly matched layers (PML) (cf. [5]), nonreflecting boundary condition methods (cf. 

ITS]), and among others (cf. [551 133 )■ An essential ingredient for the latter approach is 
to truncate an unbounded domain to a bounded domain by imposing an exact or approximate 
nonreflecting (absorbing or transparent) boundary condition at the outer artificial boundary, 
where the NRBC is designed to prevent spurious wave reflection from the artificial boundary 
(cf. the review papers |15l 116) and the references therein). The frequency-domain approaches 
for e.g., the time-harmonic Hclmholtz problems have been intensively investigated, while the 
time-domain simulations, which are capable of capturing wide-band signals and modeling more 
general material inhomogeneities and nonlinearities (cf. [311]), have been relatively less studied. 

Although some types of NRBCs based on different principles have been proposed (see, e.g., 
[121 HI [39l |34j [35J [T71 [T5J HI]), a longstanding issue of time-domain computation is the efficient 
treatment for NRBCs that can scale and integrate well with the solver for the underlying trun- 
cated problem (cf. [38l[T9]). In practice, if an accurate NRBC is imposed, the artificial boundary 
could be placed as close as possible to the scatter that can significantly reduce the computational 
cost. In this paper, we restrict our attention to the exact NRBC on the circular or spheri- 
cal artificial boundary. One major difficulty lies in that such a NRBC is global in space and 
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time in nature, as it involves the Fourier/spherical harmonic expressions in space, and history 
dependence in time induced by a convolution. The convolution kernel, termed as nonreflecting 
boundary kernel (NRBK) in [5] , is the inverse Laplace transform of an expression that includes the 
logarithmic derivative of a modified Bessel function. The rapid computation of the NRBK and 
convolution is of independent interest. Alpert, Greengard and Hagstrom [2 proposed a rational 
approximation of the logarithmic derivative with a least square implementation, which allows for 
a reduction of the summation of the poles from 0(v) to 0(logi/log |) (where ^> 1 is the order 
of the modified Bessel function and e is a given tolerance), and a recursive convolution. Jiang 
and Greengard [301 IS] further considered some interesting applications to Schrodinger equations 
in one and two dimensions. Li |23| introduced a more accurate low order approximation of the 
three-dimensional NRBK at a slightly expensive cost, where the observation that the Laplace 
transform of the three-dimensional NRBK is exactly a rational function lies at the heart of this 
algorithm. However, in many cases, the expressions are not rational functions. For instance, the 
two-dimensional NRBK also contains the contributions from the brach-cut along the negative 
real axis (see Theorem 2.1 below). Lubich and Schadle [24] developed some fast algorithm for 
the temporal convolution with 0(N t \ogN t ) operations (over N t successive time steps) arising 
from NRBCs with non-rational expressions for other equations (e.g., Schrodinger equations and 
damped wave equations). 

In this paper, we derive an analytic formula for the NRBK based on a direct inversion of 
the Laplace transform by the residue theorem (see Theorem 2.1 below). In fact, Sofronov |35j 
presented some formulas of similar type by working on much more complicated expressions of 
the kernel in terms of Tricomi's confluent hypergeometric functions. We show that with these 
formulas, we can evaluate the temporal convolution recursively and rapidly with 0(N t ) operations 
and almost without extra memory for the history dependence. Moreover, the analytic expression 
provides a useful apparatus for the stability and convergence analysis. It is worthwhile to remark 
that Chen [5] reformulated the two-dimensional wave problem into a first-order system in time 
and showed the well-posedness of the truncated problem with an alternative formulation of the 
NRBC. 

It is known that the nonlocality of the NRBC in space can be efficiently handled by Fourier/spherical 
harmonic expansions when the scatter is a disk or a ball. Recently, a systematic approach, 
based on the boundary perturbation technique (also called the transformed field expansion (TFE) 
method (cf. [H])), has been developed in [221 [JHl EO] for time-harmonic Helmholtz equations in 
exterior domains with general bounded obstacles, under which the whole algorithm boils down 
to solving a sequence of Helmholtz equations in a 2-D annulus or a 3-D spherical shell. In 
this paper, we highlight that this notion can be extended to time-domain computation, though 
it has not been investigated before as far as we know. In this paper, we propose an efficient 
spectral-Galerkin method with Newmark's time integration for the truncated wave equations in 
an annulus or a spherical shell, and provide ample numerical results to show the efficiency of the 
solver and high accuracy of NRBC from several angles. 

The rest of the paper is organized as follows. In Section 2, we present the formulation of NR- 
BCs, and derive the analytic formulas for NRBKs. In Section 3, we present some properties of 
NRBK and analyze well-posedness of the truncated wave equation. In Section 4, we outline the 
notion of the TFE method and propose an efficient spectral-Galerkin and Newmark's time inte- 
gration scheme for the truncated wave problem in regular domains. We provide ample numerical 
results in Section 5. 
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2. Evaluation of nonreflecting boundary kernels 



In this paper, we consider the time-domain acoustic scattering problem with sound-soft bound 
ary conditions on the bounded obstacle: 

d$U = c 2 AU + F, in fi^ := R d \ D, t > 0, d = 2,3; 
U = U , d t U=U l , in Q x , t = 0; 



U = G, on r D , t>0; d t U + cd n U = o(\x 



(l-d)/2N 



oo, t > 0. 



(2.1) 
(2.2) 
(2.3) 



Here, D is a bounded obstacle (scatter) with Lipschitz boundary To, c > is a given constant, 



and the radiation condition (2.3), where n = x/\x\, corresponds to the well-known Sommerfeld 
radiation condition in the frequency domain. Assume that the data F, Uq and U± are compactly 
supported in a 2-D disk or a 3-D ball B of radius b. 

A common way is to reduce this exterior problem to the problem in a bounded domain by 
imposing an exact or approximate NRBC at the artificial boundary Tf, := dB. In what follows, 
we shall focus on the wave equation truncated by the exact circular or spherical NRBC: 

d?U = c 2 AU + F, mCl:=B\D, t>0, d = 2,3; (2.4) 

U = U , d t U = Ui, infl, t = 0; U = G, on T D , t > 0; (2.5) 

d r U = T d (U), at r = b, t > 0, (2.6) 

where Td{U) is the so-called time-domain DtN map. 



2.1. Formulation of Td{U). We first present the expression of Td{U) in (2.6), and refer to 
e.g., [13 E] (and the original references therein) for the detailed derivation. It is known that 
the problem dO-dO, exterior to D = B with F = U = V x = and G = U\ r=b (i.e., the 



Dirichlet data taken from the interior problem ( 2.4 )-( 2.6 )), can be solved analytically by using 



Laplace transform in time and separation of variables in space in polar coordinate (r, <fi) /spherical 
coordinate (r, 9, (f). By imposing the continuity of directional derivative with respect to r across 
the artificial boundary r = b, we obtain the boundary condition (2.6) with 



T d {U) 



where 



ldU 


u 


elk 


~ 2r 


1BU 


U 


c~dt 


r 




:=£ 



r—b 



r—b 



- <J n {t)*U n {b,t)e in t 
|n|=0 

oo n 
n=0 |m|=0 



/2 {i)*U nm {b,t)Y™{e,<i>), d=3, 



1 

26 



s K' v (sb/c) 



1/2. 



(2.7) 



(2- 



cK v {sb/c)_ 

Here, K v is the modified Bessel function of the second kind of order v (see, e.g., [USD]), an d 
£ _1 [_ff(s)] is the inverse Laplace transform of a Laplace transformable function h(t) with 



H(s)=£{h(t)}( s ) 



e- st h{t)dt, seC, Re(s) > 0. 



In (2.7), {Y^ 71 } are the spherical harmonics, which arc orthonormal as defined in and 
{U n }/{U nm } are the Fourier/spherical harmonic expansion coefficients of U\ r= b- Note that the 
convolution is defined as usual: (/ * g)(t) = J Q f(t — T)g{r)dT. 

Alternatively, we can represent Td(U) by the temporal convolution in terms of the expansion 
coefficients of d t U\ r =b- More precisely, we define 



uj v (t) := bJ v {t; d) :- 



(d-l)c I 1 



(2.9) 
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and note that u>' u (t) = ca v {i). Then, we find from (2.7 1 and integration by parts that 



T d (U) 



ldU 



J2 uj n (t)*d t U n (b,t)e in t 



|n|=0 

oo n 

E w «+i/2W * d t U nm (b,t)Y™(9,0), 

{n=0 \m\=0 



(2.10) 



where for d = 2,3, 

UJ u (t) = CT 1 



1 - 



(d-2)c K(sb/c) 



(t), i/ = n,n + l/2. 



(2.11) 



2bs K u (sb/c 

Hereafter, we term cr„ and ui v as the nonrefiecting boundary kernels (NRBKs). Since K^ n (z) = 
K n (z) (see Formula 9.6.6 in [1]), it suffices to consider w„ and a n with n > 0, for d = 2. 

Remark 2.1. In the expressions of o~ v andio v , some terms are added, e.g., s/c and 1/(26) in 
(2.8 1, for the purpose of removing the singular part from the ratio K' v jK v . Indeed, recall the 
asymptotic formula for fixed v > and large \z\ (see Formula 9.7.2 of pQ): 



K v {z) 



^L e -. {l+ ^ +0(2 -, )} , 



for |argz| < 3tt/2, and the recurrence relation: 

zK' v {z) = vK v {z)-zK v+1 {z). 

One verifies that 



-l - 



i 



0{z~ 2 ). □ 



(2.12) 

(2.13) 
(2.14) 



K v (z) ' 2z 

Remark 2.2. We find that the use of the NRBK u;„ is more convenient, if one reformulates 
( 2.4 ) into a first- order (with respect to the time variable) system (cf. [B]), and it is more suitable 
for analysis as well, while the NRBK o~ v is more appropriate for computation. □ 

We see that NRBCs are global in both time and space. To solve the truncated problem 
( 2.4 )-( 2.6 ) efficiently, we need to (i) invert Laplace transform to compute NRBKs; (ii) deal with 
temporal convolutions efficiently; and (iii) handle the nonlocality of the NRBC in space effectively. 
The rest of the paper will address these issues. 

2.2. Evaluation of the NRBKs. Our starting point is to invert the Laplace transform via 
evaluating the Bromwich's contour integral: 



0-v(t) = 



1 



27ri 



7+ooi 



s K v (sb/c) 
/ 7 -ooi — cK v (sb/c) 

for v = n,n + 1/2 with n > 0, where 



1 

2b 



e t8 ds = 



c 

26 2 7Ti 



7+ooi 



7— 001 



1 K'Jz) 



F v (z)e czt/b dz, (2.15) 



(2.16) 



and 7 is the Laplace convergence abscissa, which is a generic constant greater than the real part 
of any singularity of F v (z). 

In order to use the residue theorem to evaluate (2.151, we need to understand the behavior of 
the poles of F u (z), i.e., the zeros of K v (z). 

Lemma 2.1. Let v > be a real number. 

(i) If z is a zero of K v (z), then its complex conjugate z is also a zero. Moreover, all complex 
conjugate pairs of zeros lie in the second and third quadrants with Re(z) < 0. 
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(ii) The total number of zeros of K v (z) is the even integer nearest to v — 1/2, if v — 1/2 is 
not an integer, or exactly v — 1/2, if v — 1/2 is an integer. 

(iii) All zeros of K n (z) and K n+ i/ 2 {z) are simple, and lie approximately along the left half of 
the boundary of an eye-shaped domain around z = (see Figure 2.1). 



Proof. The properties (i) and (ii) can be found from Page 511 of [40] . We now consider the 
property (iii). As a consequence of (i), it suffices to consider the zeros of K v (z) in the third 
quadrant and on the negative real axis (i.e., with — 7r < argz < —tt/2) of the complex plane. 
According to Formula 9.6.4 of [T], we have the following relation between K v {z) and the Hankel 
function of the first kind: 



K v {z) 



-7r < argz < 



(2.17) 



yeWffW^), _ 2 . 

which implies that all zeros of K v (z) in the third quadrant (i.e., with — 7r < argz < —n/2) are 
obtained by rotating all zeros of Hp\z) in the fourth quadrant (i.e., with — 7r/2 < argz < 0) by 
an angle —tt/2. Recall that the zeros of Hl(z) in the fourth quadrant lie approximately along 
the boundary of an eye-shaped domain around z — (see Figure 9.6 and Page 441 of pQ), whose 
boundary curve intersects the real axis at z ~ n and the imaginary axis at z = —ina, where 
a = y/tl - 1 fa 0.66274 and to « 1.19968 is the positive root of cothi = t. □ 

For clarity, let M„ be the total number of zeros of K v (z) with v = n, n + 1/2, that is, 

{the largest even integer nearest to n — 1/2, for K n (z), 
n, for K n+ y 2 {z). 

We plot in Figure ! 



(2.18) 



2.1 



some samples of zeros of K n {z) (and K n+ i/ 2 (z) for various n, and visualize 
that for a given n, the zeros sit on the left half boundary of an eye-shaped domain that intersects 
the imaginary axis approximately at ±ni, and the negative real axis at —na with a ps 0.66274 
(see the dashed coordinate grids) as predicted by Lemma 2.1 (iii). 





FIGURE 2.1. Distributions of {zj}j% of K n (left) and K n+1 / 2 (right) for various 



With the above understanding of the poles of the integrand F v (z) in (2.16), we now present 
the exact formula for the NRBKs a v {t) with v = n, n+ 1/2. 

Theorem 2.1. Let v = n, n + 1/2 with n > 0, and let {zjj^f^ be the zeros of K v (z). Then 
• for d = 2, 

-ctr/b 



o- n (t) 



c 



Kl(r)+^P n {r) 



dr 



(2.19) 
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• for d = 3, 



where I n (z) is the modified Bessel function of the first kind (cf. pQ). 



(2.20) 



We sketch the proof in Appendix|A]by applying the residue theorem to the Bromwich's contour 
integral (2.151. We remark that Sofronov |35] derived formulas of similar type by working on 
much more complicated expressions in terms of Tricomi's confluent hypergeometric functions. 
However, the formulas in the above theorem are more compact and informative. 

Remark 2.3. Based on a delicate study of the logarithmic derivative of the Hankel function 
Hi 1] (z), Alpert et al. g] ( see Theorem 4.1 and Lemma 4.2 in [5]) obtained the following formula: 



1 



iz 

2 



1 N v 



3 = 1 



7r cos(i/7r) 



(2.21) 



-dr, 



7ri J cos 2 [v-K)K 2 {r) + [nl v (r) + sm(ini)K l/ (r)) 2 ir + z 

for any v ^ n+ 1/2, where h u \, ftv,2j •• • , n u,N v are zeros of Hy(z), which number N v . Interest- 
ingly. (2.191 can be derived from (2.21), which is justified in Appendix |B) □ 



Remark 2.4. Thanks to (2.91, we obtain from Theorem 2.1 the expression of L) v (t) 
• for d= 2, 

Mn irm -, —ctr/b 



for d = 3, 



c 
2b 



V (e ctz ^ b - 1) + (-1)" / 

3=1 J ° 



1 



e 



r{K*(r)+**I*(r)} 



dr 



,ctz"/b _ 



l), i/ = n + l/2. □ 



(2.22) 



(2.23) 



2.3. Computation of the improper integral in (|2.19|). The computation of the two-dimensional 

(2.24) 



NRBK requires to evaluate the improper integral involving the kernel function 

1 1 



W n (r) 



n > 0, r > 0, 



whose important properties are characterized below. 

Lemma 2.2. For any n > and any real r > 0, we have 

(i) G n (r) is a convex function of r, and W n (r) attains its maximum at a unique point 

(ii) For large n, we have the uniform asymptotic estimate: 

n 



W n {nn) 



-sech(2n9) := W n (nn), 



for k > 0, where 



9 = 6(/c) := a/1 + K 2 + In - 



(2.25) 



(2.26) 



i + Vi + 

Approximately, the maximum value of W n (r) attains at r — na with a sa 0.66274 being 
the root of O, and the maximum value is approximately nyl + a 2 /it « 0.38187n. 
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Proof, (i). We find from Page 374 of [T] that for a given n, K n (r),I n (r) > 0, and K n (r) 
(resp. In{r)) is monotonically descending (resp. ascending) with respect to r. From the series 
representation (see Formula 9.6.10 of [I]): 

J2k 



E 

k=0 



2 2k k\{n + k)V 



n > 0, 



we conclude that I%(r) > 0. Moreover, since K n (r) satisfies (see Formula 9.6.1 of [T]) 

r 2 Kl(r) + rK' n (r) - (r 2 + n 2 )K n (r) = 0, 
we have K'^(r) > 0. Therefore, a direct calculation shows that G' r [(r) > 0, so G n (r) is convex. 



One verifies that G„(0+) = G„(+oo) = +oo for all n, which follows from (2.12), and the 
asymptotic properties (see [T] again): 



K u (r) , 
for < r -c 1, and 



In r, 
2 V2 



if y = 0, 
if v > 0, 



v > 0, 



1 

27rr 



r > 1, f > 0. 



(2.27) 



(2.28) 



Since G n (r) is convex, G n {r) attains its minimum at a unique point tq. Thanks to G' n {r) 
~ W n (r) /W 2 (r), W n (r) has a unique maximum at the same point tq. 
(ii). Recall that for large n (see Formulas (9.7.7)-(9.7.11) of [J): 



K n {nn) 



-nQ 



n0 



2n (I + k 2 ) 1 /*- 



I n (nit) 



I' n (nn) 



/^(I + k 2 ) 1 / 4 ' 

1 (l + K 2 ) 1 / 4 
. 1 

/27rn K 



(2.29) 



which, together with (2.24), leads to the asymptotic estimate (2.25). Thus, the maximum value 
of W n (r) approximately attains at the unique root of 0(re), which turns out to be a sa 0.66274 
as in Figure 2.1 and the maximum value is about ny/l + a 2 /it ps 0.38187n. □ 



W n (r) 

Asymptotics 

(na, n(1 +a 2 ) 1/2 /it) 



20 25 30 35 



FIGURE 2.2. Left: graph of ©(«),« e (0,8] defined in (2.26 1. Right: W„ against the 



asymptotic estimate W n ("•") for samples of n = 5, 15, 30, 45. 



We depict in Figure 2.2 (left) the graph of 0(«), and highlight the zero point (a, 0). Observe 



that 0(k) grows like «, when k > a. We also plot in Figure 2.2 (right) several sample graphs of 
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W n (solid lines) and the asymptotic estimate W n ("•") for n = 5, 15, 30, 45, and particularly mark 
the asymptotic point (no, ny/l + a 2 /tt) of the maximum of W n (r) obtained in Lemma 2.2 (ii). 



Observe that even for small n, the asymptotic estimate provides a very accurate approximation of 
W n . These properties greatly facilitate the computation of the improper integral. Indeed, we can 
truncate (0, oo) to a narrow interval (of length about 20 for n > 5 in our numerical computations 
in Section [5]) around the point r = na. 

2.4. Rapid evaluation of the temporal convolution. Remarkably, the presence of the time 
variable t in the exponentials in (2.19) and (2.20) allows us to eliminate the burden of history 
dependence of the temporal convolution easily. More precisely, given a function g(t), we define 



f(t;r) 



-ctr /b 



g(t) 



-c(t — r)r/b 



g(r)d,T. 



One verifies readily that 



f(t + At;r) =e- cAtr ' b f{t-r) + 



t+At 



-c{t+At-r)r/b 



g(r)dT, 



(2.30) 



so f(t;r) can march in t with step size At recursively. This enables us to compute the time 
convolution rapidly. For example, in the 2-D case, 



g](t) 



<7 n (i - T)g(r)dT 



(-l)"c 
b 2 

M„ 



c 

1 



Mn ft 

/ e c ^- T >V b g{ T )d7 



[jf 



-c(t~r)r/b 



3 = 1 



o Kl(r)+^P n {r)iJ 

(-l)"c '•+ 00 
b 2 



g(r)dr 



dr 



(2.31) 



/(t;r)W„(r)dr. 



Thanks to (2.30), [a n * g](t + Ai) can be computed recursively from the previous step and the 
history dependence is then narrowed down to [i, t + At]. We also refer to Subsection 4.4 for more 
detailed discussions. 



3. A PRIORI ESTIMATES 



In this section, we analyze the well-posedness of the truncated problem (2.4)-(2.6), and provide 
a priori estimates for its solution. 

We first recall the Plancherel or Parseval results for the Laplace transform. 

Lemma 3.1. Let s = Si +1S2 with si,S2 S K- If fig are Laplace transformable, then 



1 

2^ 



C[f}(s)C[g}( S )ds 2 



-2s x t 



f(t)g(t)dt, V Sl > 7 , 



(3.1) 



where 7 is the absissa of convergence for both f and g, and g is the complex conjugate of g. 
Proof. This identity can be proved by following the same lines as for (2.46) in 8 (or [TT]). □ 
For notational convenience, we introduce the modified spherical Bessel function (cf. [40 ): 



k n{z) = \j — K n+X / 2 {Z) 



and by (2.11), 



w n +i/a(t) = £ 1 



k„(z) 
k n (sb/c) 



1 

~2~z 



K, 



+1/2 



(i), n>0. 



(3.2) 
(3.3) 
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The following properties are also indispensable for the analysis. 

Lemma 3.2. Let s = si + is 2 with s±, s 2 G Then for any s\ > 0, 

Z' n {sb/c) 
Z n (sb/c) 



(i) Re(, 



ZnWc) 
Z n (sb/c) 



< 0, Re( 



Z n (sb/c) 



<0; 



(ii) Im^. 



< 0, Vs 2 > 0, 



(3.4) 



where Z n (z) = K n (z) or k n (z). 



Proof. The results with Z n = K n were proved in Chen [BJ. We next prove (3.4) with Z n = k n 
by using a similar argument. By applying Laplace transform to (2.1)-(2.3), exterior to D = B 

(3.5) 



with F = Uq = Ui = and G = U\ r -b, and denoting u — C[U], we obtain 

- c 2 Au + s 2 u = 0, in il cxt = R d \B, s £ C, Re(s) > 
u\ r = b = ij; = £[U\ r=b ]; cd r u + su = o(r {1 - d)/2 ), d=2,3, 

which admits the series solution 

k n (sr/c) 



i(r, 0, </», s) = 2^ 



fc„ (sb/c, 

n=0 nK 1 ' |m|=0 



(3.6) 



where {f„ m } are the spherical harmonic expansion coefficients of tp. Multiplying the first equa- 
tion of (3.5) by u and integrating over f2b jP := B p \ B, where B p is a ball of radius p > b, the 
imaginary part of the resulting equation reads 



2sis 2 



\u\ 2 dx - c 2 Im 



du 

■m b . dn 



ud'y = 0, 



where n is the unit outer normal of dtt 



b, P - 



Since Si > 0, multiplying (3.7) by s 2 yields 

du 



S2 Im 



{r=h} 



Or 



ud'y < s 2 Im 



{r=p} 



du 
Or 



u d~f. 



(3.7) 



(3.8) 



It is clear that u = k n (sr/c)/k n (sb/c) , ip nm (s)Y^ ri (6,4>) is a solution of (3.5), so using the ortho 

k' n {sp/c) 



onality of {i^ 1 }, we obtain from (3.8) that 
b 



-Im s 2 s 



KWc) 

k n (sb/c) 



\1> 



i:m 2 < -Im(s 2 s 
c V 



k n (sb/c) 



IVv, 



(3.9) 



We find from (2.12) and (3.2) that k' n (sp/c) decays exponentially if s± > 0, so the right hand 

30. Thus 

Imfs 2 s- 



side of (3.8) tends to zero as p —> +00. Thus, letting p —> +00 in (3.9) leads to 

k' n {sb/ c y 



k n {sb/c) 



< 0. 



(3.10) 



If S2 > 0, we obtain (ii) in (3.4). Next, we prove the first inequality of (i) in (3.4). Recall the 



formula (see Lemma 2.3 of [H]): 

,2 = 1 

2 



\K n+1/2 (sr)f 



-Hog 



l 2 ^ 2 „ 2 _l^2 



dr 



K n+x /2{r) — , si > 0, 



which implies that |if n+1 / 2 (sr)| 2 is monotonically descending with respect to r. The property 
f r \K n+1/2 {sr)\ 2 < 0, together with @, implies Rc(sg^ 
7i + i72 with 71,72 £ R, we know from (3.10) that 71 < and s 2 7 2 < 0. Therefore, 

71 + 172 \ 1 



< 0. Denoting fl ^^/cj 



Re 



k n (sb/c) 



r(si7i + s 2 7 2 ) < 0. 



This ends the proof. 



□ 
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With the above preparations, we can derive the following important property. 
Theorem 3.1. For any v <E L 2 (0,T), we have 

w v * v](t)v(t)dt < { \v(t)\ 2 dt, VT>0, n>0, 



(3.11) 



where uj v is the NRBK given by (2.11) (or (3.3) ford = 3). 



Proof. Let v = vl, 0T] , where 1 [0 T] is the characteristic function of [0, T]. Then we obtain from 



( [XT] ) that for d = 3, 

r-T 



II 



u n+ i/2 * v](t)v(t)dt 



+ 00 



Ht [uJ n+1/2 *v]{t)i{t)dt 



1 

1 

2n 



+°° rW 



K(sb/c) 
k n (sb/c) 



\C[v](s)\ ds2 



k n (sb/c) 



C[v](s)\ ds 2 



2tt 



|£[u](s)| ds 2 - 



It is clear that by (3.1 1 with / = g = v, 



i r+ co „ r+°° 



v(t)\ 2 dt = 



v{t)\ 2 dt. 



Using (2.13) and the properties K v (z) = K v {z) and £[w](s) = £[u](s), we find 



1 

2^ 



K(sb/c) 
k n (sb/c) 



C[v](s)\ ds2 



Re 



fk!Jjsb/c)\, |2 



k n (sb/c 



Re 



/ k' n (sb/c) \ : 
\kJsblc)) 11 



C 



(3.12) 



v(r)dT (s) 



<k n (sb/c) 

Thus, we conclude from Lemma 13.21 and the above identities that for Si > 0, 



dsq 



,-2si* 



Hi+1/2 * 1#M*)cft < 



|u(t)| 2 dt. 



(3.13) 



Notice that the asymptotic formulas in (2.27) are also valid for complex r (see Formula 9.6.9 of 
PQ), which, together with Lemma |2~Tj implies that k' n (sb/c)/k n (sb/c) is analytic for all Re(s) > 
but \s\ 7^ 0, and lim Sl _j. + \s\ 2 k' n (sb/c)/k n (sb/c) exists for all s 2 > 0. Hence, the integral in (3.12) 
is finite. By letting si — > + in (3.13) leads to the desired result for d = 3. 

The result (3.11 ) with d — 2 can be proved in a similar fashion. 

Corollary 3.1. Suppose that v' 6 L 2 (0,T) luzf/i i>(0) = 0. TTien we /iai>e 



□ 



1 



\v'{t)\ 2 dt- 



d-1 
46 



«(r)| s 



VT > 0, 



(3.14) 



/or v = n, n + 1/2, where a v (t) is the NRBK defined in (2.8). 



Proof. By (2.9), we have uj' v {t) — co- u {t),ui v {Q) = — — 2 ^ c . Thus, we obtain from integration by 
parts and the fact v(0) =0 that 

(d-l)c 



'/](*) 
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By Theorem 3.1 with v' in place of v, 

'■' (d-l)c '■' 



[lj v * v'](t)v' \t)dt 
(d-l)c. 



2/; 



v(t)v' \t)dt + c [<j„ * v](t)v' (t)dt 



-lb 



\v{T)\ 2 +c [a u *v}(t)v'(t)dt< \v\t)\ 2 dt 



This gives (3.14). 



□ 



Now, we are ready to analyze the stability of the solution of the truncated problem (2.4)- 



(2.6 1. To this end, we assume that the scatter D is a simply connected domain with Lipschitz 
boundary T D , and the Dirichlet data G = on Y D . Denote X := {U £ H 1 ^) : U\ Td = 0}, and 
let (•) ')l 2 (si) an d II • Hl 2 (si) be the inner product and norm of L 2 (Vl), respectively. 



Theorem 3.2. Let U(e X) be the solution of yLty-yLty with G = 0. If U Q £ H 1 ^),^ £ 
L 2 (fl) and F £ L 1 (0, T; L 2 (fl)) for any T > 0, then we have VU £ L°°(0,T; (H 1 (SY)) d ), d t U £ 
L°°(0,T;£ 2 (O)), and there holds 

\\dtU\\L=°(V,T;L 2 (n)) + c||Vt/|| L oo( 0: T;L2(o)) 

< C(||^i||^(n) + c||Vf/ ||^(a) + \\F\\ L i {Q . T] mn))), 

where C is a positive constant independent of any functions and c, b. Moreover, if the source term 
F = 0, we have the conservation of the energy: E'(t) = for all t > 0, where 



(3.15) 



Bit) 



c 2 \VU\ 2 )dx-2c 2 



o Jr 



T d (U)d T Ud-fdT. 



(3.16) 



Proof. Multiplying (2.4 1 by 2d t U and integrating over f2, we derive from the Green's formula 
that for any t > 0, 

d 



dt 



(\\dtU\ 



i 2 (a) + c2 ||VE/|||2(n) 



2c 2 / T d (U)d t Udj = 2(F,d t U) LHQy (3.17) 
Jr* 



Integrating the above equation over (0, t), we find that for any t > 0, 



ll^lli^)+c 2 ||W||i 2(n) -2 C 2 



T d (U)d T Ud-fdi 



o 



(F, d r U) L2 (n) dr + 1 1 Ui 1 1 1 2 (n) + c 2 1 1 VC/ 1 1 h (O) 

- ) ||9t[/|| L ~ ( o,T;L 2 (n))l|F|| i i (0 ,T;L 2 (a)) + ||l/l||ia (n) + C 2 1| Vf/ 1| ^(o) • 



(3.18) 



3 (0,T;L 2 (fi)) 

We next show that for any t > 0, 



/ / T d {U)d T Ud<ydT < 0. 
Jo Jr,, 



(3.19) 



For d — 3, it follows from (2.10), Theorem 3.1 and the orthogonality of {Y^™} that 



o Jr 



T d {U)d T Ud 1 dr = -- \\d T U\\ 2 L , {Tb) dT 



oo n 



n=0 \m\=0 



~E E / [v n +i/2 *d T U nm (b,T)]d T U nm (b,T)dT 



j3~fTJ 



„t oo n „t 

\ \\dTU\\ 2 La , rb) dT+-J2 E / |5rC>„ m (6,r)| 2 dr = 0. 

JO C , m—l „J0 



n=0 | m |=0 
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This verifies (3.19) with d = 3. Similarly, one can justify (3.191 with d = 2. Consequently, the 



estimate (3.151 follows from (3.18), (3.19) and the Cauchy-Schwarz inequality. 



Letting F = in (3.18) leads to the conservation of energy. 



□ 



Remark 3.1. Note that (i) ifTd(U) is given by (2.7), we can use Corollary 3.1 to verify (3.19) 7 
and (ii) a similar estimate for d — 2 was derived by [5] with a slightly different argument. □ 



4. Spectral-Galerkin approximation and Newmark's time integration 



This section is devoted to numerical approximation of the truncated problem (2.4)-(2.6) with 



a focus on the treatment for the NRBC (2.6). Note that if the scatter D is a disk/ball (i.e., 



in (2.4)-(2.6) is an annulus/spherical shell), the nonlocality of the NRBC in space becomes 



local in the space of Fourier/spherical harmonic coefficients in polar/spherical coordinates. Cor- 
respondingly, the truncated problem can be reduced to a sequence of one-dimensional problems 
with mixed boundary conditions at the outer artificial boundary. In order to deal with a general 
scatter, one may resort to the transformed field expansion method (TFE) (cf. [H]), which has 
been successfully applied to time-harmonic Helmholtz equations (cf. [351 HH1 ISO])- The use of 



this method allows us to solve a sequence of truncated problems ( 2.4 )-( 2.6 1 in regular domains. 
Accordingly, it is essential to construct a fast and accurate solver for (2.4)-(2.6) in an annulus or 



a spherical shell, which will be the concern of this section. We shall report the combination of 
the solver with the TFE method in a future work as the implementation is rather involved. 



4.1. Notion of TFE. We outline the TFE method for (|2.4[)-(|2.6[) with d = 2 and 
{r = b + n(4>) :O<0<2^}, 



r. 



b > bo + max 

0G[O,2tt] 



where r(4>) = bo + rj(<p) is the parametric equation of the boundary of the scatter D. 
• Make a change of variables 

, (b-b )r - bq{4>) , 
" = (b-bo)-^) ' 
which maps fi to the annulus f^o = {(r',<//) : bo < r' < b, < 4>' < 27r}. To simplify 
the notation, we still use U, F, r, 6 etc. to denote the transformed functions or variables. 



Then the problem (2.4)-(2.6) becomes 

dfU = c 2 AU + F + J(ry, 17), in Q , t > 0; 
U = U , d t U = U u in n 0> t = 0; 

U\ r=bo =G, t>0; (d r U~T d (U))\ r=b = L(r,,U)\ r=b , t > 0, 

where J{n, U) and L(rj 1 U) contain differential operators with nonconstant coefficients. 
• To solve the transformed problem efficiently, we adopt the boundary perturbation tech- 
nique by viewing the obstacle as a perturbation of the disk. More precisely, we write 
n = e( and expand the solution U as U(r,<j),t;e) — Ui{r, (f), t)e l , and likewise for 

the data. Formally, we obtain the sequence of equations after collecting the terms of e l 
(see [13): 

d 2 t Ui = c 2 AU t + Ft + J(ji, ■ ■ • , Ui-i), in fio, t > 0; 
Ui = U ,i, d t Ui = U u , in n 0) f = 0; 

U l \ r=bo =G l , <>0; (d r Ui-T d {U l ))\ r=b = L(r ] ,U l _ 1 )\ r=b , t > 0. 



Solve the above equation for I = 0, 1, 
imation. 



and sum up the series by using a Pade approx- 
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4.2. Prototype equation and dimension reduction. Under this notion, the whole algorithm 
boils down to solving a sequence of prototype equations. More precisely, we consider 



d 2 t U = c 2 AU + F, in fl = {x <= R d : b < \x\ < b}, t> 0; 

U = U , d t U = U u inflo, t = 0; U\ r=bo = 0, (d r U - T d {U)) \ = 0, t> 0, 



(4.20) 



where Td(U) is the DtN map as before. We expand the solution and given data in Fourier 



series/spherical harmonic series. Then the problem (4.20), after a polar (in 2-D) and spherical 
(in 3-D) transform, reduces to a sequence of one-dimensional problems (for brevity, we use u to 
denote the Fourier/spherical harmonic expansion coefficients of U, and likewise, we use Uo,U\ 
and / to denote the expansion coefficients of Uq,Ui and F, respectively): 



d 2 u 



2 9 ^ d " 1 ?)+^n4 = /' b <r<b, i>0; 
or j r z 



v 



u\ t =o = u 



_1 dr 
du 

dt 



= u\, b <r<b; u\ r=bo — 0, t > 0; 



(4.21) 



'1 du du 
.c dt dr 



H — u] — I a u (t — r)u(b,T)dT, t > 0, 

2r / r=b In 



where /3 n — n 2 , n(n + 1) and v = n,n + 1/2 for d — 2, 3, respectively. 

Since er„ is real, the real and imaginary parts of u and the given data can be decoupled. In 
what follows, we assume they are real. 

4.3. Spectral-Galerkin approximation in space. We apply the Legendre spectral Galerkin 



method to approximate (4.21) in space. For convenience of implementation, we transform the 
interval (bo,b) to the reference interval 7 = (—1,1) by r = ^^x + ^™ with x € 7, and denote 
the transformed functions by v(x, t) = u{r,t) 1 h{x,t) — f(r,t) and Vi(x) = u t (r) with i = 0,1, 

:ai 

\(x + c o) d ~ 1 g -) + c 2 /3 n ^ - = h, xel, t>0; 

(4.22) 



respectively. Then (4.21) can be reformulated as 
d 2 v 



c" 2 



dv 



dt 2 

v(x,0) = v (x), ^(x,0)=vi(x), xel; u(-l,t)=0, i>0 

d — 1 \ r* . 

> o 



/Idv 2 dv d - 1 \ , , /"* , N , N , 



where the constants c = and cq = P^ 3 - . 



b-bn 



Let Vat := {ip € Pn ■ VK - 1) = 0}: where Pjy is the set of all algebraic polynomials of degree 
at most N. The semi-discretization Legendre spectral-Galerkin approximation of (4.22) is to find 
vn(x, t) G Vat for all f > such that for all w £ Vn, 



(vn, w)w + c(l + c ) d 1 v N (l,t)w(l) + c 2 (d x v Nl d x w) m + c 2 /3 n (v N (x + c ) 2 ,w) 

-i;jv(l,t) - <7„(f) * 1^(1,*)) io(l) = (J N h,w)„, 



2c 2 , , rf i 
(l + co)^ 1 ' 



(4.23) 



wjv(ar, 0) = uo,jv( x )> 0) = ui,jv(z), x e 7 

where (•, -) ro is the (weighted) inner product of 7^(7) with the weight function vj = (x + cq) 
denotes dfv or as usual, In is the interpolation operator on (N + l) Legendre-Gauss-Lobatto 
points, and Vq n, «i,iv € -Pat are suitable approximations of the initial values. 



d-l 



Like Theorem 3.2 we have the following a prior estimates. 
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Theorem 4.1. Let vn be the solution of (4.23|. Then for all t > 0, 

\\ d tVN\\ 2 L 2^ {I) + C 2 (]\d x V N \\\i^ + Pn\\vN/{x + Co)!)^^)) 

< c(ll u i,A r lll^(/) +c 2 {^d x VQ, N \\ 2 L 2 ij{I) + f3 n \\v ,N / (x + c )\\ 2 L ^ {I} ) + \\I N h\\ 2 L ^j^ , 
where C is a positive constant independent of N, c and b. 



(4.24) 



This estimate holds for (4.22) with v, Vq, V\, h in place of Vpf,Vo n,Vi n, Ijyh, respectively. 



Proof. Taking w = d t v^ in (4.23), and integrating the resulted equation with respect to t, we 



use Corollary |3.1| and the argument similar to that for Theorem |3.2| to derive the estimates. □ 
Remark 4.1. Equipped with Theorem\A.l\ we can analyze the convergence of the semi-discretized 



scheme (4.23) as on Page 341 of |32J. 



□ 



We next examine the linear system of (4.23). As shown in |3TJ [32], it is advantageous to 
construct basis functions satisfying the underlying homogeneous Dirichlet boundary conditions. 
Let L[(x) be the Legendre polynomial of degree I (see, e.g., [35]), and define 4>k(x) = L k (x) + 
L k+1 (x). Then (j) k (-l) = and V N = span{0 fc : < k < N - l}. Note that <f> k (l) = 2. Setting 

N-l 

vjsr(x,t) = y] Vj(t)(j)j(x), v(t) = (%,"■ ,vn-iY, rriij = (0j,<&) w , Sj 

3=0 

we obtain the system: 

N-l 



Mv + aEv + c 2 (S + f3 n M)v + fiEv - ca (a u * y~] v. 



l=h. 



(4.25) 



3=0 



with v(0) = Vq and v(0) = V\, where M — (mij),S — (sij),M = (rhij) and E = 11* is 
an N X N matrix of all ones. In (4.25), Vo and V\ are column- N vectors of the expansion 
coefficients of i>o,jv and vi.n in terms of {</>&}, and the constants a — 8c(l + co) d_1 /(& — bo) and 
t i = ic 2 (d-l)(l + c ) d - 1 '/(b(b-b )). 

4.4. Newmark's time integration. To discretize the second-order ordinary differential system 
(4.25), we resort to the implicit second-order Newmark's scheme, which has wide applications in 
the field of structural mechanics (see [27] HI]). To this end, let At be the time-stepping size, and 
let {v m ,v m ,v m } be the approximation of {v,v,i>} at t = t m = mAt, and h m — h(t m ). 
We first take care of the convolution term \a v * J^J^q 1 (*m+i) m (4.25): 



N-l N-l 

<*v * E %] (Wi) = E 

3=0 3=0 



N-l 



o-n{t m +i - T)vj{r)dT + ^ / cr n (t m+1 - T)vj(r)d 



3=0 



At 
At 



N-l 



N-l 



M°)E^ m+1 +M Ai )E^ 



3=0 

a u (0)Ev m+1 



3=0 



N-l 

E 

3=0 



o- n (t m+ i - T)vj{r)dT 



where we used the Trapezoidal rule (of order 0(At 3 )) to approximate the integral over (t m ,t m+ i), 
and denoted by g m the approximation of the remaining terms. Note that g m depends on the 
history v l (0 < I < m), but fortunately, it can be evaluated recursively and rapidly as described 
in Subsection 2.4 Moreover, only the history data ^2ij(Q < I < m) need to be stored, and no 
any matrix-vector multiplication is involved. As a result, the burden of the history dependence 
can be eliminated. 
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Denoting 



A = M, B = aE, C = c 



c 2 (S + f3 n M) + (/x-ca^a„(0))£, (4.26) 



and / m+1 = h m+1 + cag m , we carry out the full scheme for (4.221 as follows. 

(i) Set v° = Vq and v° = "0\, and compute v° from the system ( 4.25[ ) by 

v° = M-^h - aEv° ~c 2 (S + f3 n M)v° - pEv }, (4.27) 

where the convolution term vanishes at t = 0. 

(ii) For m > 0, compute v m+1 from 

(A + VAtB + 9At 2 C)i m+1 = f n+1 -B(v m + {1- d)Atv m ) 

/ l-26» \ ( 4 - 2 8) 

- C[y m + Atv m + -^—At 2 v m ^ 



1 Ofl 

v rn+l = y m + A ^m + ___ At ^ + ff^m+l^ 



and update u m+1 and u m+1 by 
• } m+1 = v m + A; 

» m+1 = v m + (1 - d)Atv m + $Atv rn+ \ (4.30) 

respectively. 

Remark 4.2. In general, the Newmark's scheme is of second- order and unconditionally stable, 
if the parameters satisfy i9 > \ and 9 > |(| + {see, e.g., |41j). 

T/ie computational cost of the full algorithm lies in solving ( 4.28[ ) with (4.26). The matrices 
M,M and S are sparse with small bandwidth under the basis {4>k}i an d their entries can be 
evaluated exactly by using the properties of Legendre polynomials. The matrices B and C in 
(4.26) appear to be full, but since all entries of E are one, the matrix-vector multiplication can 
be carried out in O(N) operations. □ 

5. Numerical results 

In this section, we present various numerical results to show the accuracy of the NRBC and 
and convergence of the proposed algorithm. 

5.1. Testing problem and setup. We examine the NRBC and the scheme via the following 
problem with an exact solution. 



Proposition 5.1. The exterior problem ( 2.1 )-( 2.3 1 with d = 2, D = {x e R 2 : |x| < b },F = 
Uq = Ux = 0, and U\ r —t, = G, admits the solution: 

oo 

U{r,<j>,t) = t4M)e in0 , r > b , [0,2tt), * > 0, (5.1) 

where 

U(rf) = l°' t</?0 ' (5 2) 

ny,J \H n (r,t)*G n (t-M + ^b /rG n (t-p ), t>(3 , 

with /3q = (r — bo)/c, G n (t) being the Fourier expansion coefficient of G(<f>,t), and 

, ^ K n+1 (rzf/b ) 
n n {r,i) > K'(z n ) 

°i=i n{ i } (5.3) 
i i i\n c f°° In(rp/bo)K n (p) — K n (rp/bo)I n (p) 

— ctpfbo j 

b J K^(p) + 7T 2 I 2 (p) 
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Proof. We sketch its derivation in Appendix [C] 

Remark 5.1. Notice that U n satisfies the one- dimension problem: 



d 2 U n c 2 d / dU n 

o ( r 

dt 2 r dr\ dr 



U n \t=o 



dt 

ldU n , dU n 



+ c 2 n 2 ^=0, b <r<b, t > 0; 



0, b < r < b; U n \ r=b(l = G n , t > 0; 



1 



a n (t-T)U n (b,T)dT, t>0. □ 



. c dt dr ' 2r 
We take the Dirichlet data: 

G(<j),t) = A ie - t((bocos *- I '» !+ ''» sin +-v»)*) sm p (ut), 



□ 

(5.4) 
(5.5) 
(5.6) 

(5.7) 



where A\ = 10, L = 0.1, x s = y s = 2.1, b = 2, and uj,p are to be specified later. We compute the 
Fourier coefficient 

g n (f) = ^sin>t) 



g-t((bo cos <j>—x 3 ) +{b sin 4>~Ve) ) e —in4 



via the fast Fourier transform with sufficient Fourier points so that the error of numerical inte- 
gration can be neglected. 

The first important issue is to compute the exact solution and o~ n very accurately. We point 



out that the convolution in (5.2) can be computed exactly, and the improper integral in (5.3) can 



be treated by a similar manner as the improper integral of o~ n {t) in Subsection 2.3 We adopt a 
Newton's iteration method to compute the zeros of K n (z) accurately for moderate large n, and 
use a (composite) Legendre-Gauss quadrature to evaluate the integral with high precision. 



For the readers' reference, we tabulate in Table 5.1 some samples of a n (t) with b = 3, c = 5 



TABLE 5.1. Some samples of a„(t). 



n 


<Jn(0.1) 


<J„(2) 


n 


(Tn(O.l) 


<Tn(2) 





5.922023678764810e-2 


1.127355852971518c-2 


5 


-5.394296605508512 


2.652114252130534e-3 


1 


-1.770775252065292e-l 


-1.849758830245009c-2 


6 


-7.505798337812085 


-1.628640292879232c-3 


2 


-8.766785435408012c-l 


-3.516167498630298c-3 


7 


-9.786660490985828 


1.032611050422504e-3 


3 


-2.012014067208546 


6.186828109004313e-3 


8 


-12.14241492386835 


-6.793719076424341e-4 


4 


-3.538153424430011 


-4.319820513209325e-3 


9 


-14.47330651171318 


4.628600071971645e-4 



We plot in Figure [5TT] the exact solution with p = 2, bo = 2, b = 4, cj = 107r at various t. We 
see that the Dirichlet data G(<p, t) acts as a dynamical wave-maker. 



> o 




FIGURE 5.1. Time evolution of exact solution with oj = 107T. 
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5.2. Accuracy of the NRBC. Now, we examine accuracy of the NRBC. Notice that 

• 1 dU n dU, 



max \d r U M ~ T d (U M ) \ < 



M 

E • 

|n|=0 



c dt 



dr 



Un 

2r 



(T n (t)*U n (b,t) 



where T d (U M ) is defined in (J^Yf with d = 2, and U M (b,(f>,t) = J2\n\=oUn{b,t)e in<t > is the trun- 
cation of the exact solut ion i n Proposition |5.1| We choose the same parameters as those for the 
exact solution in Figure 5.1 and take M — 32 so that {\U n (b, t)\} are sufficiently small for all 



modes \n\ < M and all t of interest. Note that the differentiations in e n (t) can be performed 
exactly on the exact solution. 

In Table 15.21 we tabulate the errors: 



max e n (t), 

\n\<M 



E [ M\t) 



M 

E 

|n|=0 



e n (t), 



at some samples of t, compute the errors at the outer artificial boundary r — b, where the exact 
solution has the magnitude as large as possible. 



TABLE 5.2. The errors E$(t) and E$(t). 



t 


u = 10ir,b = 2.22 


oj = 107r,6 = 2.75 


w = 20tt, 6 = 2.38 


w = 20tt, b = 2.87 


















0.5 


2.637c-16 


4.902c-16 


7.026c-17 


1.148e-16 


6. 245c- 17 


1.603c-16 


8.327c-17 


1.752c-16 


1 


4.684c-16 


8.406c-16 


2.099c-16 


2.630c-16 


9.021e-17 


2.390c-16 


2. 272c-16 


4.405c-16 


5 


2.567c-16 


4.386c-16 


1.131e-15 


1.176o-15 


7.251c-16 


9.053c-16 


1.318o-15 


1.473C-15 


10 


7.772c-16 


1.020c- 15 


2.207c-15 


2.227o-15 


1.457c-15 


1.960o-15 


2.387c-15 


2.575e-15 



We see from Table 5.2 that in all tests, the errors are extremely small. This validates the 



formula for a n in Theorem 2.1 and high accuracy of the numerical treatment. 



5.3. Numerical tests for the spectral-Galerkin-Newmark's scheme. Hereafter, we pro- 
vide some numerical results to illustrate the convergence of the numerical scheme described 
in Section |4j In the following computations, we take the Dirichlet data given by (5.7) with 
A\ = 10, l = 0.1, x s = y s = 2.1 and b = 2 as before. 

Given a cut-off number M > 0, we compute the numerical solutions {U^} of ( 5.4 1-( 5.6 ) for 
the modes |n| < M by using the spectral-Galerkin and Newmark's time integration scheme. 
The full numerical solution of the problem in Proposition 
(r,t)e m ^ ', and the exact solution Um(t, 4 1 , t) = 



5.1 



is then defined as Uj^{r,<f>,t) 



M=o' 



2| n |_ U n (r,i)e in 4> is evaluated as 



before. Once again, we choose M such that the Fourier coefficient of the exact solution \U n \ < 
1Q-16 f or | n | y m. The numerical errors are measured by 

\U%{;t)-U n {;t)\ vN M - „„„ II^AT, 



it) 



max 
|n|<Af 



\L 2 ,N> 



(t) 



max \\U"(;t) 

\n\<M " 



,N> 



where || • ||l 2 .jv is the discrete i 2 -norm associated with the Legendre-Gauss-Lobatto interpolation, 
and || • ||i«>,iv is the corresponding discrete maximum norm. 

To test the second-order convergence of the Newmark's scheme, we choose N — 50 so that the 



error of the spatial discretization is negligible. We provide in Table 5.3 the numerical errors and 
the order of convergence for various t with M — 15,6 = b,u = n, and p = 6. As expected, we 
observe a second-order convergence of the time integration. 

Next, we fix the time step size At = 10 -5 and choose different N to check the accuracy 
in spatial discretization. The convergence behavior is illustrated in Table 5.3 With a small 
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TABLE 5.3. The errors E$[(t) and E^(t), and order of convergence. 



t 


At 




order 




order 


t 


At 




order 




order 




lc-03 


1.21e-05 




1.93e-05 






lc-03 


1.23c-05 




2.01c-05 




1 


5e-04 


3.02e-06 


2.00005 


4.83e-06 


1.99994 


3 


5e-04 


3.07e-06 


2.00007 


5.02c-06 


1.99997 




le-04 


1.21c-07 


1.99993 


1.93e-07 


1.99986 




le-04 


1.23e-07 


1.99997 


2.01c-07 


1.99998 




5e-05 


3.02e-08 


2.00012 


4.83e-08 


1.99906 




5e-05 


3.07e-08 


2.00007 


5.02c-08 


2.00004 




lc-03 


1.23e-05 




2.01e-05 






lc-03 


1.23c-05 




2.01c-05 




2 


5e-04 


3.07e-06 


1.99998 


5.02e-06 


1.99997 


4 


5e-04 


3.07e-06 


2.00007 


5.02e-06 


1.99997 




le-04 


1.23e-07 


2.00001 


2.01c-07 


1.99998 




le-04 


1.23c-07 


1.99997 


2.01c-07 


1.99998 




5e-05 


3.07e-08 


2.00002 


5.02e-08 


1.99993 




5e-05 


3.07e-08 


2.00007 


5.02c-08 


2.00004 



t=2.5 




t=4 



b=5 




(a) 



(b) 



(c) 



FIGURE 5.2. Propagation of exact and numerical solutions. 

number of modes in space, we observe a fast decay of the errors, which is typical for the spectral 
approximation. 



TABLE 5.4. The errors E$(t) and E$(t). 



t 


N = 8 


N = 10 


AT = 16 


N = 32 


Em 






E^it) 


E^(t) 


Em 


Em 


Em 


0.5 


2.349c-04 


2.687c-04 


2.253e-05 


2.846e-05 


8.224c-07 


1.207c-06 


9.229c-09 


1.639e-08 


1.0 


3.877c-04 


4.179c-04 


1.708e-05 


2.042e-05 


1.833c-07 


2.777c-07 


1.959c-09 


2.673c-09 


1.5 


3.276c-04 


3.596e-04 


5.965e-06 


7.146e-06 


3.072c-08 


5.238c-08 


1.508e-09 


1.762e-09 


2.0 


4.105c-04 


4.368e-04 


1.092e-05 


1.188c-05 


7.353c-09 


1.238e-08 


1.237c-09 


2.119c-09 


2.5 


3.239c-04 


3.607e-04 


5.724e-06 


6.997c-06 


2.381c-09 


3.182c-09 


1.488c-09 


1.734c-09 


3.0 


4.113e-04 


4.380e-04 


1.095e-05 


1.189e-05 


1.404e-09 


2.893c-09 


1.227c-09 


2.009c-09 


3.5 


3.238e-04 


3.608e-04 


5.747e-06 


7.036e-06 


1.554e-09 


2.063c-09 


1.488e-09 


1.732e-09 


4.0 


4.113e-04 


4.380e-04 


1.095c-05 


1.189e-05 


1.269c-09 


2.061c-09 


1.227c-09 


2.007c-09 



Finally, we plot in Figure 5.2 the numerical solution (red nets) against the exact solution (blue 
smooth surface), where the exact solution in (a) and (b) is in the annulus: 2 < r < 10, while 
the numerical solution within 2 < r < 5. Figure [5~2{ c) shows the wave propagation through the 
artificial boundary at b = 5, < <f> < 2w for < t < 4 (numerical solution) against < t < 4.5 
(exact solution). We see that the red nets and blue surfaces agree well in these plots. Moreover, 
we observe the waves pass the boundary transparently. This shows that the proposed scheme is 
very stable and has high-order accuracy. 



Concluding remarks 
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We proposed in this paper analytic and accurate numerical means for the time-domain wave 
propagation with exact and global nonreflecting boundary conditions. We derived the analytic 
expressions of the convolution kernels and presented highly accurate methods for their evalu- 
ations. We analyzed the stability of the truncated problem and provided efficient numerical 
schemes. Ample numerical results were given to demonstrate the features of the method. We 
shall report the combination of the methods with the boundary perturbation technique in a fu- 
ture work. The techniques and ideas in this paper will be useful to study the Maxwell equations 
and elastic wave propagations. 



Appendix A. Proof of Theorem 12.1 



We first consider d — 2. Note that F n (z) is a multi- valued function with branch points at 
z = and infinity. The contour L for (2.15) is depicted in Figure A.l| (left) with the branch-cut 
along the negative real axis. 





FIGURE A.l. Contour L for the inverse Laplace transform. Left: d — 2. Right: d — 3. 



We know from Lemma 2.1 that F n (z) has a finite number of simple poles in the second and 
third quadrants, Therefore, for any n > 0, it follows from the residue theorem that 



lim <b F n (z)e czt/b dz = 2m V Res \F n (z) 

«->+«> J L -h 

r->0 + J_i 



e czt/ \z: 



zt/b 



1 K'Jz) 
z H h z —, { 

2 K n {z) 



3=1 



e ctz ? /b z] 1 . 



Thus, by |2.15|) and 
2b 2 ni 



a n (t) = 2tti Ve ctz ? /fc z! 1 - lim [/ F n (z)e czt/b dz + [ F n {z)e czt/b dz 

Ipi fl^+oo [Jgg Jap 

3- L r-y0 + 

+ f^F n (z)e czt/b dz+ [_F n (z)e czt/b dz + [ F n {z)e czt ' b dz 
Jdef Jfg Jg2 



(A.l) 



(A.2) 



In view of (2.14), we find from the Jordan's lemma (cf. [TTJ [5] ) and a direct calculation that 

0, (A.3) 



lim 

.R-S-+00 
r->0+ 



_F n {z)e czt/b dz+ F n {z)e czt/b dz+ lF n {z)e czt/b dz 
bS Jdef Jg2 
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since each contour integral tends to zero. Thus, we only need to evaluate the integrals along the 
line segments CD and FG. We have 

~K' n (re-™) K' n (re™y 



lim 

r^0+ 



CD 



F n (z)e czt/b dz + / F n (z)e czt/b dz 



FG 



-ctr/b 



K n (re-™) K n (re™) 



dr. 



By Formula 9.6.31 of pQ, 

K n (re*) = e- n ^K n (r) - ml n {r), K n {re~^) = e n ^K n (r) + nil n (r), 



which, together with (2.13), implies 



K' n (rj*) K' n {re-™) = K n+1 (re- iv ) K n+1 (re™) 
K n {re™) K n (re-™) K n {re~™) " K n (re™) 

= eM^r) + ml n+1 (r) e^ n + l ^K n+l {r) - ^i/ w+1 (r) 
e n7ri X„(r) + 7riJ n (r) e- n7rl iT n (r) - 7ri/ n (r) 

_ (-l)*27ri[X- B+ i (r)I„(r) + tf„(r)I„ +1 (r)] _ (-l)"2yri 



(A.4) 



Kl{r)+^Il{r) r[Kl{r) + -K*Il{rj\ 

where we used the Wronskian identity (see Formula 9.6.15 of [1]): 

I n (z)K n+1 (z) + I n+1 (z)K n (z) = z^ 1 . 



A combination of (A.1)-(A.4) leads to 

M„ 



°n(t) = u v z'v* 2 ?/ 6 + / ~ ^ \ re z'l M 



(A.5) 



(A.6) 



which is the expression (2.19). 

Now, we turn to the three dimensional case. It is important to notice that the kernel function 
F n+ i/2(z) is not a multi- valued complex function, as opposed to the two dimensional case. Indeed, 
although K n+ i/2(z) is multi-valued, in view of the formula (see Page 80 of |40) ) : 

w») = v^t H ' ( "-^V - Vn£0 ' (A - 7) 

the fact \j \fz can be eliminated from the ratio K' n+l j 2 / K n+ i/2- Thanks to this property, we use 
the contour in Figure A.l (right), by the residue theorem and Jordan's lemma, we have 

/■7+ooi 



lim J> F n+l/2 {z)e czt ' b dz 

R-^+ooJ L 

M„ 



F n+l/2 (z)e czt ' b dz 



(A. 



2th ^ Res [F n+1/2 (z)e czt / b ,z»] = 2th ]T 



3=1 



for v = n + 1/2. This leads to the formula (2.20). 



□ 



Appendix B. Justification for Remark 12.3 



Thanks to (2.17), we have 



HP\z) . K(-iz) 



'K v {-\ z y 



Thus, by (2.21 ) with v = n and —\z = bs/c, we find 

N„ 



1 sK' n (bs/c) _ ic 



2b 



E 



"n,3 



K n (bs/c) V £j s + b 2 Jo K&r) + 7r»J*(r) f + 



(-1)" 



-dr. 
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Applying the inverse Laplace transform to both sides of the above identity leads to 



c- 1 



s i s K' n (bs/c) 
c 26 cK n {bs/c) 



N„ 



f_]\»»g— Ct/f> 



-dr, 



VJ Kl{r)+^P n {r) 

which turns out to be (A. 6), since — ih n> j are zeros of K n (z) and N n = M n . □ 

Appendix C. Proof of Proposition 15.11 



This problem can be solved by Laplace transform and separation of variables. Indeed, the 



Fourier coefficients {U n } in (5.1) can be expressed as 

-if K n (sr/c) 



U n (r,t)=£- 



{t)*G n (t), r>b . 



(C.l) 



<K n (sb a /c), 

We next sketch the evaluation of the inverse Laplace transform by using the residue theorem as 

g " (sr/c) ~- • F>fe- sfio , where f3 = (r - b Q )/c. In 



in Appendix 



A 



for a n {t). Using (2.12), we find 



K n (sb /c) 



order to use the Jordan's lemma (cf. [HJ[E]), we write 

\K n (sb /c)) [t) L \ E V K n (sb /c) Sri) 



c- 1 



-Pos 



(t) 



H n (r,t - p )U Po (t) 



(C.2) 



6(t-0o), 



where H n (r,t) = C^ 1 (e' 3 °' i ^"(sbijl) ~ V^o/ r ) (t), Up (t) is the unit step function (which takes 
1 for t > f3 , and for t < /3o), and S(t) is the Dirac delta function. By applying the residue 
theorem and Jordan's lemma to the Bromwich's contour integral: 



H n {r,t) 



27ri 



7+ooi 



7—001 



.0 



j s K n {sr/c) 
K n (sb Q /c) 



e ts ds, 



(C.3) 



we obtain (5.3) by using the contour in Figure A.l (left) and the same argument as in Appendix 
|a| Denoting H n (r,t) := H n (r,t — /3 ), the expression (5.2) follows from (C.2). We leave the 
details to the interested readers. □ 
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